require(stringi)
require(estimatr)
require(effectsize)
require(FindIt)
require(RMeCab)
require(quanteda)
require(textir)

factor.list <- c("group", "attribute", "sphere", "definition")
factor.labels <- c("Group", "Attribute", "Sphere and mechanism", "Definition")
level.labels <- list(c("Advantaged", "Disadvantaged"), 
                     c("Gender", "Nationality", "Sexuality", "Education"), 
                     c("Public (taste-based)", "Public (stereotypical)", 
                       "Public (statistical)", "Public (customer needs)", "Private"), 
                     c("None", "Minimal", "Detailed"))

#### ChatGPT survey data ####
GPT.data <- read.csv("GPT_4_data.csv")

## convert ChatGPT responses into numerical outcome values
GPT.data$rating <- GPT.data$rating.raw <- 
  as.numeric(stri_replace_all_regex(stri_sub(stri_trans_general(GPT.data$answer, "Fullwidth-Halfwidth"), 
                                             1, stri_locate_first_regex(GPT.data$answer, "理由")[, 1] - 1), 
                                    "\\n| |\\(.+\\)|\\:|【|】|答え", ""))

## manual correction of irregular or incomplete responses
GPT.data$answer[which(is.na(GPT.data$rating.raw))]

GPT.data$rating[which(is.na(GPT.data$rating.raw))[1]] <- NA
GPT.data$rating[which(is.na(GPT.data$rating.raw))[2]] <- 1.5
GPT.data$rating[which(is.na(GPT.data$rating.raw))[3]] <- 5
GPT.data$rating[which(is.na(GPT.data$rating.raw))[4]] <- 6
GPT.data$rating[which(is.na(GPT.data$rating.raw))[5]] <- 4
GPT.data$rating[which(is.na(GPT.data$rating.raw))[6]] <- 4
GPT.data$rating[which(is.na(GPT.data$rating.raw))[7]] <- NA
GPT.data$rating[which(is.na(GPT.data$rating.raw))[8]] <- 5.5
GPT.data$rating[which(is.na(GPT.data$rating.raw))[9]] <- 5.5
GPT.data$rating[which(is.na(GPT.data$rating.raw))[10]] <- 3.5
GPT.data$rating[which(is.na(GPT.data$rating.raw))[11]] <- 3
GPT.data$rating[which(is.na(GPT.data$rating.raw))[12]] <- 5
GPT.data$rating[which(is.na(GPT.data$rating.raw))[13]] <- 6
GPT.data$rating[which(is.na(GPT.data$rating.raw))[14]] <- 6
GPT.data$rating[which(is.na(GPT.data$rating.raw))[15]] <- 4.5
GPT.data$rating[which(is.na(GPT.data$rating.raw))[16]] <- NA
GPT.data$rating[which(is.na(GPT.data$rating.raw))[17]] <- 1.5
GPT.data$rating[which(is.na(GPT.data$rating.raw))[18]] <- 3.5
GPT.data$rating[which(is.na(GPT.data$rating.raw))[19]] <- 4

## center outcome values to compute marginal means
GPT.data$rating.centered <- GPT.data$rating - mean(GPT.data$rating, na.rm = TRUE)

## create variables corresponding to experimental conditions
GPT.data$attribute <- factor(ifelse(GPT.data$Q < 11, "gender", 
                                    ifelse(GPT.data$Q < 21, "nationality", 
                                           ifelse(GPT.data$Q < 31, "sexuality", "education"))), 
                             levels = c("education", "gender", "nationality", "sexuality"))
GPT.data$attribute.id <- ifelse(GPT.data$attribute == "education", 4, 
                                as.numeric(GPT.data$attribute) - 1)
GPT.data$attribute.reordered <- factor(GPT.data$attribute, 
                                       levels = c("gender", "nationality", "sexuality", "education"))
GPT.data$disadvantaged <- ifelse(GPT.data$attribute %in% c("gender", "education"), 
                                 ifelse(GPT.data$Q %% 2 == 1, 1, 0), 
                                 ifelse(GPT.data$Q %% 2 == 0, 1, 0))
GPT.data$advantaged <- 1 - GPT.data$disadvantaged
GPT.data$group <- factor(ifelse(GPT.data$disadvantaged == 1, 
                                "disadvantaged", "advantaged"), 
                         levels = c("advantaged", "disadvantaged"))
GPT.data$sphere <- factor(ifelse(GPT.data$Q %% 10 == 1 | GPT.data$Q %% 10 == 2, "taste", 
                                 ifelse(GPT.data$Q %% 10 == 3 | GPT.data$Q %% 10 == 4, "stereotype", 
                                        ifelse(GPT.data$Q %% 10 == 5 | GPT.data$Q %% 10 == 6, "statistics", 
                                               ifelse(GPT.data$Q %% 10 == 7 | GPT.data$Q %% 10 == 8, "needs", "private")))), 
                          levels = c("private", "taste", "stereotype", "statistics", "needs"))
GPT.data$sphere.id <- ifelse(GPT.data$sphere == "private", 5, 
                             as.numeric(GPT.data$sphere) - 1)
GPT.data$sphere.stereotype.ref <- relevel(GPT.data$sphere, "stereotype")
GPT.data$sphere.statistics.ref <- relevel(GPT.data$sphere, "statistics")
GPT.data$sphere.needs.ref <- relevel(GPT.data$sphere, "needs")
GPT.data$sphere.reordered <- factor(GPT.data$sphere, 
                                    levels = c("taste", "stereotype", "statistics", "needs", "private"))
GPT.data$gender <- (GPT.data$attribute == "gender") * 1
GPT.data$nationality <- (GPT.data$attribute == "nationality") * 1
GPT.data$sexuality <- (GPT.data$attribute == "sexuality") * 1
GPT.data$education <- (GPT.data$attribute == "education") * 1
GPT.data$taste <- (GPT.data$sphere == "taste") * 1
GPT.data$stereotype <- (GPT.data$sphere == "stereotype") * 1
GPT.data$statistics <- (GPT.data$sphere == "statistics") * 1
GPT.data$needs <- (GPT.data$sphere == "needs") * 1
GPT.data$private <- (GPT.data$sphere == "private") * 1
GPT.data$definition <- factor(c("none", "minimal", "detailed")[GPT.data$C], 
                              levels = c("none", "minimal", "detailed"))
GPT.data$none <- (GPT.data$definition == "none") * 1
GPT.data$minimal <- (GPT.data$definition == "minimal") * 1
GPT.data$detailed <- (GPT.data$definition == "detailed") * 1
GPT.data$definition.ref <- relevel(GPT.data$definition, "minimal")

#### human survey data ####
human.raw.data <- read.csv("human_data.csv")

## reshape the data to long format
human.data <- data.frame(id = rep(1:nrow(human.raw.data), 4), 
                         thread_id = rep(1:4, each = nrow(human.raw.data)), 
                         Q = c(human.raw.data$Q.1, human.raw.data$Q.2, 
                               human.raw.data$Q.3, human.raw.data$Q.4), 
                         rating = c(human.raw.data$rating.1, human.raw.data$rating.2, 
                                    human.raw.data$rating.3, human.raw.data$rating.4), 
                         text = c(human.raw.data$text.1, human.raw.data$text.2, 
                                  human.raw.data$text.3, human.raw.data$text.4), 
                         C = rep(human.raw.data$C, 4))

## centered outcome values to compute centered marginal means
human.data$rating.centered <- human.data$rating - mean(human.data$rating)

## create variables corresponding to experimental conditions
human.data$attribute <- factor(ifelse(human.data$Q < 11, "gender", 
                                      ifelse(human.data$Q < 21, "nationality", 
                                             ifelse(human.data$Q < 31, "sexuality", "education"))), 
                               levels = c("education", "gender", "nationality", "sexuality"))
human.data$attribute.id <- ifelse(human.data$attribute == "education", 4, 
                                  as.numeric(human.data$attribute) - 1)
human.data$attribute.reordered <- factor(human.data$attribute, 
                                         levels = c("gender", "nationality", "sexuality", "education"))
human.data$disadvantaged <- ifelse(human.data$attribute %in% c("gender", "education"), 
                                   ifelse(human.data$Q %% 2 == 1, 1, 0), 
                                   ifelse(human.data$Q %% 2 == 0, 1, 0))
human.data$advantaged <- 1 - human.data$disadvantaged
human.data$group <- factor(ifelse(human.data$disadvantaged == 1, 
                                  "disadvantaged", "advantaged"), 
                           levels = c("advantaged", "disadvantaged"))
human.data$sphere <- factor(ifelse(human.data$Q %% 10 == 1 | human.data$Q %% 10 == 2, "taste", 
                                   ifelse(human.data$Q %% 10 == 3 | human.data$Q %% 10 == 4, "stereotype", 
                                          ifelse(human.data$Q %% 10 == 5 | human.data$Q %% 10 == 6, "statistics", 
                                                 ifelse(human.data$Q %% 10 == 7 | human.data$Q %% 10 == 8, "needs", "private")))), 
                            levels = c("private", "taste", "stereotype", "statistics", "needs"))
human.data$sphere.id <- ifelse(human.data$sphere == "private", 5, 
                               as.numeric(human.data$sphere) - 1)
human.data$sphere.stereotype.ref <- relevel(human.data$sphere, "stereotype")
human.data$sphere.statistics.ref <- relevel(human.data$sphere, "statistics")
human.data$sphere.needs.ref <- relevel(human.data$sphere, "needs")
human.data$sphere.reordered <- factor(human.data$sphere, 
                                      levels = c("taste", "stereotype", "statistics", "needs", "private"))
human.data$gender <- (human.data$attribute == "gender") * 1
human.data$nationality <- (human.data$attribute == "nationality") * 1
human.data$sexuality <- (human.data$attribute == "sexuality") * 1
human.data$education <- (human.data$attribute == "education") * 1
human.data$taste <- (human.data$sphere == "taste") * 1
human.data$stereotype <- (human.data$sphere == "stereotype") * 1
human.data$statistics <- (human.data$sphere == "statistics") * 1
human.data$needs <- (human.data$sphere == "needs") * 1
human.data$private <- (human.data$sphere == "private") * 1
human.data$definition <- factor(c("none", "minimal", "detailed")[human.data$C], 
                               levels = c("none", "minimal", "detailed"))
human.data$none <- (human.data$definition == "none") * 1
human.data$minimal <- (human.data$definition == "minimal") * 1
human.data$detailed <- (human.data$definition == "detailed") * 1
human.data$definition.ref <- relevel(human.data$definition, "minimal")

#### distribution of demographic variables ####
## load and preprocess 2020 census data
census.data <- read.csv("census_data.csv", na.strings = "-")

census.data$sex <- c("men", "women")[census.data$cat01_code]
census.data$age <- census.data$cat02_code + 14
census.data <- subset(census.data, 18 <= age & age <= 69)
census.data$cohort <- ifelse(census.data$age < 29, "20s", 
                             ifelse(census.data$age < 39, "30s", 
                                    ifelse(census.data$age < 49, "40s", 
                                           ifelse(census.data$age < 59, "50s", 
                                                  ifelse(census.data$age < 69, "60s", "70s")))))
census.data$edu <- ifelse(census.data$cat03_code == 4 | census.data$cat03_code == 17, NA, 
                          ifelse(census.data$cat03_code == 3 | census.data$cat03_code == 11 | census.data$cat03_code == 12 | census.data$cat03_code == 13, "low", 
                                 ifelse(census.data$cat03_code == 14, "middle", "high")))
census.data$edu <- factor(census.data$edu, 
                          levels = c("low", "middle", "high"))
census.data$region <- ifelse(census.data$area_code < 8000, "Hokkaido/Tohoku", 
                             ifelse(census.data$area_code < 15000, "Kanto", 
                                    ifelse(census.data$area_code < 25000, "Chubu", 
                                           ifelse(census.data$area_code < 31000, "Kinki", 
                                                  ifelse(census.data$area_code < 40000, "Chugoku/Shikoku", "Kyushu")))))
census.data$region <- factor(census.data$region, 
                             levels = c("Hokkaido/Tohoku", "Kanto", "Chubu", "Kinki", "Chugoku/Shikoku", "Kyushu"))

census.data$value.2 <- census.data$value
census.data$prop <- census.data$value / sum(census.data$value, na.rm = TRUE)

census.data$value.2[is.na(census.data$edu)] <- NA
census.data$prop.2 <- census.data$value.2 / sum(census.data$value.2, na.rm = TRUE)

## gender distribution
# survey
round(prop.table(table(human.raw.data$gender)), 3)
# census
round(sum(census.data$prop[census.data$sex == "men"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$sex == "women"], na.rm = TRUE), 3)

## age distribution
# survey
round(prop.table(c(sum(human.raw.data$age < 30), 
                   sum(human.raw.data$age >= 30 & human.raw.data$age < 40), 
                   sum(human.raw.data$age >= 40 & human.raw.data$age < 50), 
                   sum(human.raw.data$age >= 50 & human.raw.data$age < 60), 
                   sum(human.raw.data$age >= 60))), 3)
# census
round(sum(census.data$prop[census.data$cohort == "20s"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$cohort == "30s"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$cohort == "40s"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$cohort == "50s"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$cohort == "60s"], na.rm = TRUE), 3)

## educational attainment distribution
# survey
round(prop.table(c(sum(human.raw.data$edu < 3), 
                   sum(human.raw.data$edu >= 3 & human.raw.data$edu < 6), 
                   sum(human.raw.data$edu >= 6))), 3)
# census
round(sum(census.data$prop.2[census.data$edu == "low"], na.rm = TRUE), 3)
round(sum(census.data$prop.2[census.data$edu == "middle"], na.rm = TRUE), 3)
round(sum(census.data$prop.2[census.data$edu == "high"], na.rm = TRUE), 3)

# difference in the proportion of respondents who graduated high school or less
round(prop.table(c(sum(human.raw.data$edu < 3), 
                   sum(human.raw.data$edu >= 3 & human.raw.data$edu < 6), 
                   sum(human.raw.data$edu >= 6)))[1] - 
        sum(census.data$prop.2[census.data$edu == "low"], na.rm = TRUE), 2)

## region of residence
# survey
round(prop.table(c(sum(human.raw.data$pref %in% c(2, 3, 8, 12, 16, 24, 45)), 
                   sum(human.raw.data$pref %in% c(4, 10, 14, 19, 35, 39, 41)), 
                   sum(human.raw.data$pref %in% c(1, 6, 9, 15, 26, 29, 38, 43, 47)), 
                   sum(human.raw.data$pref %in% c(13, 22, 23, 28, 33, 36, 44)), 
                   sum(human.raw.data$pref %in% c(5, 11, 17, 20, 31, 37, 40, 42, 46)), 
                   sum(human.raw.data$pref %in% c(7, 18, 21, 25, 27, 30, 32, 34)))), 3)
# census
round(sum(census.data$prop[census.data$region == "Hokkaido/Tohoku"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$region == "Kanto"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$region == "Chubu"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$region == "Kinki"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$region == "Chugoku/Shikoku"], na.rm = TRUE), 3)
round(sum(census.data$prop[census.data$region == "Kyushu"], na.rm = TRUE), 3)

#### distribution of survey ratings ####
# function to compute marginal means
MM.function <- function(data) {
  lm.result <- lm_robust(formula(rating ~ 1), data, clusters = id)
  MM.out <- cbind(mean = lm.result$coefficients, 
                  lower = lm.result$conf.low, 
                  upper = lm.result$conf.high)
  rownames(MM.out) <- names(lm.result$coefficients)
  MM.out
}

# function to compute distributions
dist.function <- function(x) {
  prop.table(table(factor(round(x), levels = 1:6)))
}

## marginal means in the ChatGPT survey
plot.MM.GPT <- matrix(NA, 17, 3)
plot.MM.GPT[1, ] <- MM.function(GPT.data)
for (i in 1:4) {
  plot.MM.GPT[i + 1, ] <- MM.function(subset(GPT.data, attribute.id == i & advantaged == 1))
}
for (i in 1:4) {
  plot.MM.GPT[i + 5, ] <- MM.function(subset(GPT.data, attribute.id == i & disadvantaged == 1))
}
for (i in 1:5) {
  plot.MM.GPT[i + 9, ] <- MM.function(subset(GPT.data, sphere.id == i))
}
for (i in 1:3) {
  plot.MM.GPT[i + 14, ] <- MM.function(subset(GPT.data, C == i))
}

## distribution in the ChatGPT survey
plot.dist.GPT <- matrix(NA, 17, 6)
plot.dist.GPT[1, ] <- dist.function(GPT.data$rating)
for (i in 1:4) {
  plot.dist.GPT[i + 1, ] <- dist.function(subset(GPT.data, attribute.id == i & advantaged == 1)$rating)
}
for (i in 1:4) {
  plot.dist.GPT[i + 5, ] <- dist.function(subset(GPT.data, attribute.id == i & disadvantaged == 1)$rating)
}
for (i in 1:5) {
  plot.dist.GPT[i + 9, ] <- dist.function(subset(GPT.data, sphere.id == i)$rating)
}
for (i in 1:3) {
  plot.dist.GPT[i + 14, ] <- dist.function(subset(GPT.data, C == i)$rating)
}

round(plot.dist.GPT, 2)

# percentage of responses on the 'not discriminatory' side to the customer needs scenario
round(sum(plot.dist.GPT[13, 1:3]), 2)

## marginal means in the human survey
plot.MM.human <- matrix(NA, 17, 3)
plot.MM.human[1, ] <- MM.function(human.data)
for (i in 1:4) {
  plot.MM.human[i + 1, ] <- MM.function(subset(human.data, attribute.id == i & advantaged == 1))
}
for (i in 1:4) {
  plot.MM.human[i + 5, ] <- MM.function(subset(human.data, attribute.id == i & disadvantaged == 1))
}
for (i in 1:5) {
  plot.MM.human[i + 9, ] <- MM.function(subset(human.data, sphere.id == i))
}
for (i in 1:3) {
  plot.MM.human[i + 14, ] <- MM.function(subset(human.data, C == i))
}

## distribution in the human survey
plot.dist.human <- matrix(NA, 17, 6)
plot.dist.human[1, ] <- dist.function(human.data$rating)
for (i in 1:4) {
  plot.dist.human[i + 1, ] <- dist.function(subset(human.data, attribute.id == i & advantaged == 1)$rating)
}
for (i in 1:4) {
  plot.dist.human[i + 5, ] <- dist.function(subset(human.data, attribute.id == i & disadvantaged == 1)$rating)
}
for (i in 1:5) {
  plot.dist.human[i + 9, ] <- dist.function(subset(human.data, sphere.id == i)$rating)
}
for (i in 1:3) {
  plot.dist.human[i + 14, ] <- dist.function(subset(human.data, C == i)$rating)
}

round(plot.dist.human, 2)

# percentage of responses on the 'not discriminatory' side to the customer needs scenario
round(sum(plot.dist.human[13, 1:3]), 2)

## Figure 1
main.labels <- c("Overall", "Male", "Japanese", "Heterosexuals", "UTokyo grads", 
                 "Female", "Chinese", "Homosexuals", "NEPU grads", 
                 "Public\n(taste-based)", "Public\n(stereotypical)", 
                 "Public\n(statistical)", "Public\n(customer needs)", "Private", 
                 "None", "Minimal", "Detailed")
png("Figure_1.png", 
    width = 6, height = 6, unit = "in", pointsize = 8, res = 2000)
layout(matrix(c(1, 18:20, 2:14, 21:23, 15:17, 24), 6, 4, byrow = TRUE))
par(mar = c(1.5, 0.5, 2, 0.5), omi = c(0, 1, 0, 0), lwd = 0.5, xpd = TRUE)
for (i in 1:17) {
  plot(NULL, NULL, bty = "n", xlim = c(0.5, 6.5), ylim = c(-1, 1), 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  segments(0.25, 0, 6.75, 0, col = "darkgray")
  segments(0.25, c(-0.8, -0.6, -0.4, -0.2, 0.2, 0.4, 0.6, 0.8), 
           6.75, c(-0.8, -0.6, -0.4, -0.2, 0.2, 0.4, 0.6, 0.8), lty = 3, col = "lightgray")
  segments(1:6, -0.95, 1:6, 0.95, lty = 3, col = "lightgray")
  for (j in 1:6) {
    polygon(c(j - 0.5, j + 0.5, j + 0.5, j - 0.5), 
            c(0, 0, plot.dist.GPT[i, j], plot.dist.GPT[i, j]), 
            col = "gray")
  }
  segments(plot.MM.GPT[i, 2], 0.1, plot.MM.GPT[i, 3], 0.1)
  points(plot.MM.GPT[i, 1], 0.1, pch = 19, cex = 0.8)
  for (j in 1:6) {
    polygon(c(j - 0.5, j + 0.5, j + 0.5, j - 0.5), 
            c(0, 0, -1 * plot.dist.human[i, j], -1 * plot.dist.human[i, j]), 
            col = "white")
  }
  segments(plot.MM.human[i, 2], -0.1, plot.MM.human[i, 3], -0.1)
  points(plot.MM.human[i, 1], -0.1, pch = 21, bg = "white", cex = 0.8)
  mtext(main.labels[i], line = 0.5, font = 2)
  if (i == 1) {
    text(6, plot.dist.GPT[1, 6], "ChatGPT", pos = 3)
    text(6, -1 * plot.dist.human[1, 6], "Human", pos = 1)
    mtext(c("Not\ndiscriminatory", "Discriminatory"), 
          1, c(0.5, -0.5), at = c(0.5, 6.5), adj = c(0, 1), cex = 0.6)
  }
}
mtext("Overall", 2, 2, outer = TRUE, at = 11 / 12, font = 2, las = 1)
mtext("Group status\nand attribute", 2, 2, outer = TRUE, at = 0.75, font = 2, las = 1)
mtext("Sphere and\nmechanism", 2, 2, outer = TRUE, at = 5 / 12, font = 2, las = 1)
mtext("Definition", 2, 2, outer = TRUE, at = 1 / 12, font = 2, las = 1)
dev.off()

#### analysis of ChatGPT survey responses ####
# H1, H2, H4, and H5
GPT.result.base <- lm_robust(rating ~ disadvantaged + attribute + sphere + definition, 
                             GPT.data, clusters = id)
round(summary(GPT.result.base)$coefficient[, c(1, 2, 4)], 3)

print(hedges_g(rating ~ factor(disadvantaged), data = GPT.data), digits = 3)
print(hedges_g(rating ~ attribute, 
               data = subset(GPT.data, attribute.id %in% c(1, 4))), 
      digits = 3)
print(hedges_g(rating ~ attribute, 
               data = subset(GPT.data, attribute.id %in% c(2, 4))), 
      digits = 3)
print(hedges_g(rating ~ attribute, 
               data = subset(GPT.data, attribute.id %in% c(3, 4))), 
      digits = 3)
print(hedges_g(rating ~ sphere, 
               data = subset(GPT.data, sphere.id %in% c(1, 5))), 
      digits = 3)

# H3: taste-based and stereotypical
GPT.result.stereotype.ref <- lm_robust(rating ~ attribute + sphere.stereotype.ref + disadvantaged + definition, 
                                       GPT.data, clusters = id)
round(summary(GPT.result.stereotype.ref)$coefficient[, c(1, 2, 4)], 3)

# H3: stereotypical and statistics
GPT.result.statistics.ref <- lm_robust(rating ~ attribute + sphere.statistics.ref + disadvantaged + definition, 
                                       GPT.data, clusters = id)
round(summary(GPT.result.statistics.ref)$coefficient[, c(1, 2, 4)], 3)

# H3: statistics and customer needs
GPT.result.needs.ref <- lm_robust(rating ~ attribute + sphere.needs.ref + disadvantaged + definition, 
                                  GPT.data, clusters = id)
round(summary(GPT.result.needs.ref)$coefficient[, c(1, 2, 4)], 3)

# robustness check for H5 (excluding education factor)
GPT.result.wo.education <- lm_robust(rating ~ disadvantaged + attribute + sphere + definition, 
                                     GPT.data, subset = attribute.id != 4, clusters = id)
round(summary(GPT.result.wo.education)$coefficient[, c(1, 2, 4)], 3)

# H6
GPT.result.H6 <- lm_robust(rating ~ disadvantaged + attribute + definition, 
                           GPT.data, subset = statistics == 1, clusters = id)
round(summary(GPT.result.H6)$coefficient[, c(1, 2, 4)], 3)

print(hedges_g(rating ~ definition, 
               data = subset(GPT.data, statistics == 1 & C %in% c(1, 3))), 
      digits = 3)

GPT.result.H6.ref <- lm_robust(rating ~ disadvantaged + attribute + definition.ref, 
                               GPT.data, subset = statistics == 1, clusters = id)
round(summary(GPT.result.H6.ref)$coefficient[, c(1, 2, 4)], 3)

print(hedges_g(rating ~ definition, 
               data = subset(GPT.data, statistics == 1 & C %in% c(2, 3))), 
      digits = 3)

GPT.result.H6.rev <- lm_robust(rating ~ disadvantaged + attribute + definition, 
                               GPT.data, subset = statistics == 0, clusters = id)
round(summary(GPT.result.H6.rev)$coefficient[, c(1, 2, 4)], 3)

print(hedges_g(rating ~ definition, 
               data = subset(GPT.data, statistics == 0 & C %in% c(1, 3))), 
      digits = 3)

GPT.result.H6.rev.ref <- lm_robust(rating ~ disadvantaged + attribute + definition.ref, 
                                   GPT.data, subset = statistics == 0, clusters = id)
round(summary(GPT.result.H6.rev.ref)$coefficient[, c(1, 2, 4)], 3)

print(hedges_g(rating ~ definition, 
               data = subset(GPT.data, statistics == 0 & C %in% c(2, 3))), 
      digits = 3)

#### analysis of human survey responses ####
# H1, H2, H4, and H5
human.result.base <- lm_robust(rating ~ disadvantaged + attribute + sphere + definition, 
                               human.data, clusters = id)
round(summary(human.result.base)$coefficient[, c(1, 2, 4)], 3)

print(hedges_g(rating ~ disadvantaged, data = human.data), digits = 3)
print(hedges_g(rating ~ attribute, 
               data = subset(human.data, attribute.id %in% c(1, 4))), 
      digits = 3)
print(hedges_g(rating ~ attribute, 
               data = subset(human.data, attribute.id %in% c(2, 4))), 
      digits = 3)
print(hedges_g(rating ~ attribute, 
               data = subset(human.data, attribute.id %in% c(3, 4))), 
      digits = 3)
print(hedges_g(rating ~ sphere, 
               data = subset(human.data, sphere.id %in% c(1, 5))), 
      digits = 3)

# H3: taste-based and stereotypical
human.result.stereotype.ref <- lm_robust(rating ~ attribute + sphere.stereotype.ref + disadvantaged + definition, 
                                         human.data, clusters = id)
round(summary(human.result.stereotype.ref)$coefficient[, c(1, 2, 4)], 3)

# H3: stereotypical and statistics
human.result.statistics.ref <- lm_robust(rating ~ attribute + sphere.statistics.ref + disadvantaged + definition, 
                                         human.data, clusters = id)
round(summary(human.result.statistics.ref)$coefficient[, c(1, 2, 4)], 3)

# H3: statistics and customer needs
human.result.needs.ref <- lm_robust(rating ~ attribute + sphere.needs.ref + disadvantaged + definition, 
                                    human.data, clusters = id)
round(summary(human.result.needs.ref)$coefficient[, c(1, 2, 4)], 3)

# robustness check for H5 (without education)
human.result.wo.education <- lm_robust(rating ~ disadvantaged + attribute + sphere + definition, 
                                       human.data, subset = attribute.id != 4, clusters = id)
round(summary(human.result.wo.education)$coefficient[, c(1, 2, 4)], 3)

# H6
human.result.H6 <- lm_robust(rating ~ disadvantaged + attribute + definition, 
                             human.data, subset = statistics == 1, clusters = id)
round(summary(human.result.H6)$coefficient[, c(1, 2, 4)], 3)

human.result.H6.ref <- lm_robust(rating ~ disadvantaged + attribute + definition.ref, 
                                 human.data, subset = statistics == 1, clusters = id)
round(summary(human.result.H6.ref)$coefficient[, c(1, 2, 4)], 3)

human.result.H6.rev <- lm_robust(rating ~ disadvantaged + attribute + definition, 
                                 human.data, subset = statistics == 0, clusters = id)
round(summary(human.result.H6.rev)$coefficient[, c(1, 2, 4)], 3)

human.result.H6.rev.ref <- lm_robust(rating ~ disadvantaged + attribute + definition.ref, 
                                     human.data, subset = statistics == 0, clusters = id)
round(summary(human.result.H6.rev.ref)$coefficient[, c(1, 2, 4)], 3)

## Figure 2
png("Figure_2.png", 
    width = 6, height = 4, unit = "in", pointsize = 8, res = 2000)
par(mar = c(3.5, 0, 3, 0), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-3.8, 5.1), ylim = c(0.5, 16), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = c(0, 3.4), col = "darkgray")
abline(v = c(-1, -0.5, 0.5, 1, 1.5, 2.4, 2.9, 3.9, 4.4, 4.9), 
       lty = 3, col = "darkgray")
segments(-1.2, c(14.5, 12:10, seq(7.5, 4.5, -1), 2:1), 
         1.7, c(14.5, 12:10, seq(7.5, 4.5, -1), 2:1), 
         lty = 3, col = "darkgray")
segments(2.2, c(14.5, 12:10, seq(7.5, 4.5, -1), 2:1), 
         5.1, c(14.5, 12:10, seq(7.5, 4.5, -1), 2:1), 
         lty = 3, col = "darkgray")
text(-4, 16, "Group status", pos = 4, font = 2)
text(-4, 15.5, "(ref: advantaged)", pos = 4)
text(-3.8, 14.5, expression(paste("Disadvantaged (", beta, ")")), pos = 4)
segments(GPT.result.base$conf.low[2], 14.5, 
         GPT.result.base$conf.high[2], 14.5)
points(GPT.result.base$coefficients[2], 14.5, pch = 19)
segments(human.result.base$conf.low[2] + 3.4, 14.5, 
         human.result.base$conf.high[2] + 3.4, 14.5)
points(human.result.base$coefficients[2] + 3.4, 14.5, pch = 21, bg = "white")
text(-4, 13.5, "Attribute", pos = 4, font = 2)
text(-4, 13, "(ref: education)", pos = 4)
text(-3.8, 12:10, 
     c(expression(paste("Gender (", gamma[1], ")")), 
       expression(paste("Nationality (", gamma[2], ")")), 
       expression(paste("Sexuality (", gamma[3], ")"))), 
     pos = 4)
segments(GPT.result.base$conf.low[3:5], 12:10, 
         GPT.result.base$conf.high[3:5], 12:10)
points(GPT.result.base$coefficients[3:5], 12:10, pch = 19)
segments(human.result.base$conf.low[3:5] + 3.4, 12:10, 
         human.result.base$conf.high[3:5] + 3.4, 12:10)
points(human.result.base$coefficients[3:5] + 3.4, 12:10, pch = 21, bg = "white")
text(-4, 9, "Sphere and mechanism", pos = 4, font = 2)
text(-4, 8.5, "(ref: private)", pos = 4)
text(-3.8, seq(7.5, 4.5, -1), 
     c(expression(paste("Public (taste-based) (", delta[1], ")")), 
       expression(paste("Public (stereotypical) (", delta[2], ")")), 
       expression(paste("Public (statistical) (", delta[3], ")")), 
       expression(paste("Public (customer needs) (", delta[4], ")"))), 
     pos = 4)
segments(GPT.result.base$conf.low[6:9], seq(7.5, 4.5, -1), 
         GPT.result.base$conf.high[6:9], seq(7.5, 4.5, -1))
points(GPT.result.base$coefficients[6:9], seq(7.5, 4.5, -1), pch = 19)
segments(human.result.base$conf.low[6:9] + 3.4, seq(7.5, 4.5, -1), 
         human.result.base$conf.high[6:9] + 3.4, seq(7.5, 4.5, -1))
points(human.result.base$coefficients[6:9] + 3.4, seq(7.5, 4.5, -1), pch = 21, bg = "white")
text(-4, 3.5, "Definition", pos = 4, font = 2)
text(-4, 3, "(ref: no definition)", pos = 4)
text(-3.8, 2:1, 
     c(expression(paste("Minimal (", zeta[1], ")")), 
       expression(paste("Detailed (", zeta[2], ")"))), 
     pos = 4)
segments(GPT.result.base$conf.low[10:11], 2:1,  
         GPT.result.base$conf.high[10:11], 2:1)
points(GPT.result.base$coefficients[10:11], 2:1, pch = 19)
segments(human.result.base$conf.low[10:11] + 3.4, 2:1,  
         human.result.base$conf.high[10:11] + 3.4, 2:1)
points(human.result.base$coefficients[10:11] + 3.4, 2:1, pch = 21, bg = "white")
axis(1, at = seq(-1, 1.5, 0.5), 
     labels = sprintf("%.1f", seq(-1, 1.5, 0.5)), lwd = 0.5)
axis(1, at = seq(2.4, 4.9, 0.5), 
     labels = sprintf("%.1f", seq(-1, 1.5, 0.5)), lwd = 0.5)
mtext("Average marginal component effect", 1, 2.5, at = c(0.25, 3.65))
mtext(c("ChatGPT", "Human"), line = 1, at = c(0.25, 3.65), cex = 1.2, font = 2)
dev.off()

#### comparison of ChatGPT and human survey results ####
## combine ChatGPT and human survey data into a single dataset
combined.data <- rbind(cbind(id = GPT.data$id, 
                             GPT.data[, c("rating", "attribute", "sphere", 
                                          "disadvantaged", "definition")], 
                             GPT = 1), 
                       cbind(id = human.data$id + 10000, 
                             human.data[, c("rating", "attribute", "sphere", 
                                            "disadvantaged", "definition")], 
                             GPT = 0))

## regression analysis with a ChatGPT indicator variable
combined.result <- lm_robust(rating ~ disadvantaged + attribute + sphere + definition + GPT, 
                             combined.data, clusters = id)
round(summary(combined.result)$coefficients[, c(1, 2 ,4)], 3)
print(hedges_g(rating ~ factor(GPT), data = combined.data), digits = 3)

## joint-significance F-tests for interactions between factors and a ChatGPT indicator variable
anova(lm(rating ~ disadvantaged + attribute + sphere + definition + GPT, 
         data = combined.data), 
      lm(rating ~ disadvantaged * GPT + attribute + sphere + definition, 
         data = combined.data))
anova(lm(rating ~ disadvantaged + attribute + sphere + definition + GPT, 
         data = combined.data), 
      lm(rating ~ disadvantaged + attribute * GPT + sphere + definition, 
         data = combined.data))
anova(lm(rating ~ disadvantaged + attribute + sphere + definition + GPT, 
         data = combined.data), 
      lm(rating ~ disadvantaged + attribute + sphere * GPT + definition, 
         data = combined.data))
anova(lm(rating ~ disadvantaged + attribute + sphere + definition + GPT, 
         data = combined.data), 
      lm(rating ~ disadvantaged + attribute + sphere + definition * GPT, 
         data = combined.data))

## compute centered marginal means for ChatGPT and human data
# function to compute contered marginal means
centered.MM.function <- function(data) {
  lm.result <- lm_robust(formula(rating.centered ~ 1), data, clusters = id)
  MM.out <- cbind(mean = lm.result$coefficients, 
                  lower = lm.result$conf.low, 
                  upper = lm.result$conf.high)
  rownames(MM.out) <- names(lm.result$coefficients)
  MM.out
}

## centered marginal means in the ChatGPT survey
plot.centered.MM.GPT <- matrix(NA, 14, 3)
plot.centered.MM.GPT[1, ] <- centered.MM.function(subset(GPT.data, advantaged == 1))
plot.centered.MM.GPT[2, ] <- centered.MM.function(subset(GPT.data, disadvantaged == 1))
for (i in 1:4) {
  plot.centered.MM.GPT[i + 2, ] <- centered.MM.function(subset(GPT.data, attribute.id == i))
}
for (i in 1:5) {
  plot.centered.MM.GPT[i + 6, ] <- centered.MM.function(subset(GPT.data, sphere.id == i))
}
for (i in 1:3) {
  plot.centered.MM.GPT[i + 11, ] <- centered.MM.function(subset(GPT.data, C == i))
}

## centered marginal means in the human survey
plot.centered.MM.human <- matrix(NA, 14, 3)
plot.centered.MM.human[1, ] <- centered.MM.function(subset(human.data, advantaged == 1))
plot.centered.MM.human[2, ] <- centered.MM.function(subset(human.data, disadvantaged == 1))
for (i in 1:4) {
  plot.centered.MM.human[i + 2, ] <- centered.MM.function(subset(human.data, attribute.id == i))
}
for (i in 1:5) {
  plot.centered.MM.human[i + 6, ] <- centered.MM.function(subset(human.data, sphere.id == i))
}
for (i in 1:3) {
  plot.centered.MM.human[i + 11, ] <- centered.MM.function(subset(human.data, C == i))
}

## Figure 3
png("Figure_3.png", 
    width = 4, height = 4, unit = "in", pointsize = 8, res = 2000)
par(mar = c(3.5, 0, 0, 0), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(-2.5, 0.9), ylim = c(0.5, 18), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = 0, col = "darkgray")
abline(v = c(-1.2, -0.9, -0.6, -0.3, 0.3, 0.6, 0.9), 
       lty = 3, col = "darkgray")
segments(-1.25, c(17:16, 14:11, 9:5, 3:1), 
         0.95, c(17:16, 14:11, 9:5, 3:1), 
         lty = 3, col = "darkgray")
text(-2.6, 18, "Group status", pos = 4, font = 2)
text(-2.5, 17:16, level.labels[[1]], pos = 4)
segments(plot.centered.MM.GPT[1:2, 2], 17:16 + 0.15, 
         plot.centered.MM.GPT[1:2, 3], 17:16 + 0.15)
points(plot.centered.MM.GPT[1:2, 1], 17:16 + 0.15, pch = 19)
segments(plot.centered.MM.human[1:2, 2], 17:16 - 0.15, 
         plot.centered.MM.human[1:2, 3], 17:16 - 0.15)
points(plot.centered.MM.human[1:2, 1], 17:16 - 0.15, pch = 21, bg = "white")
text(-2.6, 15, "Attribute", pos = 4, font = 2)
text(-2.5, 14:11, level.labels[[2]], pos = 4)
segments(plot.centered.MM.GPT[3:6, 2], 14:11 + 0.15, 
         plot.centered.MM.GPT[3:6, 3], 14:11 + 0.15)
points(plot.centered.MM.GPT[3:6, 1], 14:11 + 0.15, pch = 19)
segments(plot.centered.MM.human[3:6, 2], 14:11 - 0.15, 
         plot.centered.MM.human[3:6, 3], 14:11 - 0.15)
points(plot.centered.MM.human[3:6, 1], 14:11 - 0.15, pch = 21, bg = "white")
text(-2.6, 10, "Sphere and mechanism", pos = 4, font = 2)
text(-2.5, 9:5, level.labels[[3]], pos = 4)
segments(plot.centered.MM.GPT[7:11, 2], 9:5 + 0.15, 
         plot.centered.MM.GPT[7:11, 3], 9:5 + 0.15)
points(plot.centered.MM.GPT[7:11, 1], 9:5 + 0.15, pch = 19)
segments(plot.centered.MM.human[7:11, 2], 9:5 - 0.15, 
         plot.centered.MM.human[7:11, 3], 9:5 - 0.15)
points(plot.centered.MM.human[7:11, 1], 9:5 - 0.15, pch = 21, bg = "white")
text(-2.6, 4, "Definition", pos = 4, font = 2)
text(-2.5, 3:1, level.labels[[4]], pos = 4)
segments(plot.centered.MM.GPT[12:14, 2], 3:1 + 0.15, 
         plot.centered.MM.GPT[12:14, 3], 3:1 + 0.15)
points(plot.centered.MM.GPT[12:14, 1], 3:1 + 0.15, pch = 19)
segments(plot.centered.MM.human[12:14, 2], 3:1 - 0.15, 
         plot.centered.MM.human[12:14, 3], 3:1 - 0.15)
points(plot.centered.MM.human[12:14, 1], 3:1 - 0.15, pch = 21, bg = "white")
par(xpd = TRUE)
legend(-2.5, par()$usr[3], legend = c("ChatGPT", "Humans"), 
       pch = c(19, 21), bg = c(NA, "white"))
axis(1, at = seq(-1.2, 0.9, 0.3), 
     labels = sprintf("%.1f", seq(-1.2, 0.9, 0.3)), lwd = 0.5)
mtext("Centered marginal mean", 1, 2.5, at = -0.3)
dev.off()

#### exploratory analysis of interaction between factors ####
# function to conduct joint-significance F-tests for interaction between factors
interaction.F.test <- function(variable.1, variable.2, data) {
  lm.model.1 <- lm(rating ~ group + attribute + sphere + definition, 
                   data = data)
  lm.model.2 <- lm(as.formula(paste0("rating ~ group + attribute + sphere + definition + ", 
                                     variable.1, " : ", variable.2)), 
                   data = data)
  anova(lm.model.1, lm.model.2)$"Pr(>F)"[2]
}

## joint-significance F-tests for interaction effects
GPT.int.F.test.p.values <- human.int.F.test.p.values <- matrix(NA, 4, 4)
rownames(GPT.int.F.test.p.values) <- colnames(GPT.int.F.test.p.values) <- 
  rownames(human.int.F.test.p.values) <- colnames(human.int.F.test.p.values) <- 
  factor.list
for (i in 1:4) {
  for (j in 2:4) {
    if (i >= j) {
      next
    }
    GPT.int.F.test.p.values[i, j] <- interaction.F.test(factor.list[i], 
                                                        factor.list[j], 
                                                        GPT.data)
    human.int.F.test.p.values[i, j] <- interaction.F.test(factor.list[i], 
                                                          factor.list[j], 
                                                          human.data)
  }
}
round(matrix(p.adjust(as.vector(GPT.int.F.test.p.values), "holm"), 4, 4), 3)
round(matrix(p.adjust(as.vector(human.int.F.test.p.values), "holm"), 4, 4), 3)

## estimate average marginal interaction effects
# remove observations with a missing outcome
GPT.data.CausalANOVA <- subset(GPT.data, ! id %in% id[which(is.na(rating))])

CausalANOVA.GPT <- CausalANOVA(rating ~ group + attribute.reordered + 
                                 sphere.reordered + definition, 
                               data = GPT.data.CausalANOVA, nway = 2, 
                               family = "gaussian", 
                               cluster = GPT.data.CausalANOVA$id)

CausalANOVA.human <- CausalANOVA(rating ~ group + attribute.reordered + 
                                   sphere.reordered + definition, 
                                 ~ group : attribute.reordered + 
                                   group : sphere.reordered + 
                                   attribute.reordered : sphere.reordered, 
                                 data = human.data, nway = 2, 
                                 family = "gaussian", 
                                 cluster = human.data$id)

## Figures S1-S6
# function to plot average marginal interaction effects
AMIE.plot <- function(variable.1, variable.2, 
                      results.GPT, order.GPT, 
                      results.human = NULL, order.human = NULL) {
  out.GPT <- results.GPT$CI.table[[order.GPT]]
  nrows <- nrow(out.GPT)
  if (is.null(results.human)) {
    plot(NULL, NULL, type = "n", bty = "n", 
         xlim = c(min(out.GPT[, 3]), max(out.GPT[, 4])), 
         ylim = c(0.5, nrows), 
         xlab = "", ylab = "", xaxt = "n", yaxt = "n")
    abline(h = nrows:1, col = "lightgray")
    abline(v = axTicks(1), lty = 3, col = "lightgray")
    segments(out.GPT[, 3], nrows:1, out.GPT[, 4], nrows:1)
    points(out.GPT[, 1], nrows:1, pch = 19)
  } else {
    out.human <- results.human$CI.table[[order.human]]
    plot(NULL, NULL, type = "n", bty = "n", 
         xlim = c(min(out.GPT[, 3], out.human[, 3]), 
                  max(out.GPT[, 4], out.human[, 4])), 
         ylim = c(0.5, nrows), 
         xlab = "", ylab = "", xaxt = "n", yaxt = "n")
    abline(h = nrows:1, col = "lightgray")
    abline(v = axTicks(1), lty = 3, col = "lightgray")
    segments(out.GPT[, 3], c(nrows:1) + 0.15, 
             out.GPT[, 4], c(nrows:1) + 0.15)
    points(out.GPT[, 1], c(nrows:1) + 0.15, pch = 19)
    segments(out.human[, 3], c(nrows:1) - 0.15, 
             out.human[, 4], c(nrows:1) - 0.15)
    points(out.human[, 1], c(nrows:1) - 0.15, 
           pch = 21, bg = "white")
  }
  axis(1, lwd = 0.5)
  mtext("Average marginal interaction effect", 1, line = 2.5)
  count <- nrows
  for (i in 1:length(level.labels[[variable.2]])) {
    for (j in 1:length(level.labels[[variable.1]])) {
      mtext(paste0(level.labels[[variable.1]][j], " \u00d7 ", 
                   level.labels[[variable.2]][i]), 
            2, 0.5, at = count, las = 2)
      count <- count - 1
    }
  }
  if (! is.null(results.human)) {
    legend(par()$usr[1] - (par()$usr[2] - par()$usr[1]) * 0.1, par()$usr[3], 
           legend = c("ChatGPT", "Humans"), 
           pch = c(19, 21), bg = c(NA, "white"), 
           xjust = 1, xpd = TRUE)
  }
}

png("Figure_S1.png", 
    width = 4.5, height = 2.5, unit = "in", pointsize = 7, res = 2000)
par(mar = c(3.5, 16.5, 0.5, 0.5), lwd = 0.5)
AMIE.plot(1, 2, CausalANOVA.GPT, 5, CausalANOVA.human, 5)
dev.off()

png("Figure_S2.png", 
    width = 4.5, height = 2.5, unit = "in", pointsize = 7, res = 2000)
par(mar = c(3.5, 16.5, 0.5, 0.5), lwd = 0.5)
AMIE.plot(1, 3, CausalANOVA.GPT, 6, CausalANOVA.human, 6)
dev.off()

png("Figure_S3.png", 
    width = 4.5, height = 2.5, unit = "in", pointsize = 7, res = 2000)
par(mar = c(3.5, 16.5, 0.5, 0.5), lwd = 0.5)
AMIE.plot(2, 3, CausalANOVA.GPT, 8, CausalANOVA.human, 7)
dev.off()

png("Figure_S4.png", 
    width = 4.5, height = 2.5, unit = "in", pointsize = 7, res = 2000)
par(mar = c(3.5, 16.5, 0.5, 0.5), lwd = 0.5)
AMIE.plot(1, 4, CausalANOVA.GPT, 7)
dev.off()

png("Figure_S5.png", 
    width = 4.5, height = 2.5, unit = "in", pointsize = 7, res = 2000)
par(mar = c(3.5, 16.5, 0.5, 0.5), lwd = 0.5)
AMIE.plot(2, 4, CausalANOVA.GPT, 9)
dev.off()

png("Figure_S6.png", 
    width = 4.5, height = 2.5, unit = "in", pointsize = 7, res = 2000)
par(mar = c(3.5, 16.5, 0.5, 0.5), lwd = 0.5)
AMIE.plot(3, 4, CausalANOVA.GPT, 10)
dev.off()

#### analysis of open-ended responses on reasons for judgments ####
## extract open-ended responses from ChatGPT outputs
GPT.data$text <- stri_replace_all_regex(stri_sub(GPT.data$answer, 
                                                 stri_locate_first_regex(GPT.data$answer, "理由")[, 2] + 1, 
                                                 stri_length(GPT.data$answer)), 
                                        "\\n|】", "")

## preprocessing of text data
text.data <- rbind(data.frame(GPT.data[, c("id", "rating", "gender", "nationality", 
                                           "sexuality", "education", "taste", 
                                           "stereotype", "statistics", "needs", "private", 
                                           "disadvantaged", "minimal", "detailed", "text")], 
                              GPT = 1), 
                   data.frame(human.data[, c("id", "rating", "gender", "nationality", 
                                             "sexuality", "education", "taste", 
                                             "stereotype", "statistics", "needs", "private", 
                                             "disadvantaged", "minimal", "detailed", "text")], 
                              GPT = 0))
text.data <- subset(text.data, ! is.na(text))

# calculate the number of characters in each response
text.data$nchar <- stri_length(text.data$text)

# tokenize text and remove functional words
text.data$mecabbed.text <- NA
for (i in 1:nrow(text.data)){
  # remove punctuations marks and symbols
  text.i <- stri_replace_all_regex(text.data$text[i], "[^ぁ-んァ-ヶーa-zA-Z0-9一-龠0-9A-Za-z]", " ")
  # unicode normalization
  text.i <- stri_trans_nfkc(text.i)
  # unify the labels of fictitious persons' names
  text.i <- stri_replace_all_regex(text.i, "Aさん|Bさん|Cさん|Dさん", "Xさん")
  # tokenize using MeCab
  mecabbed.i <- unlist(RMeCabC(text.i, 1))
  # remove functional words
  text.data$mecabbed.text[i] <- stri_c(mecabbed.i[which(! names(mecabbed.i) %in% c("助詞", "助動詞"))], collapse = " ")
}

# ChatGPT's data excluding education scenarios
GPT.mnlm.data <- subset(text.data, GPT == 1 & education == 0)
round(mean(GPT.mnlm.data$nchar))

# document-feature matrix for ChatGPT
GPT.dfm <- dfm_trim(dfm(tokens_ngrams(tokens(phrase(GPT.mnlm.data$mecabbed.text)), 1)), 
                    min_docfreq = 2)
GPT.dfm.matrix <- Matrix(as.matrix(GPT.dfm), sparse = TRUE)
dim(GPT.dfm.matrix)

# humans' data excluding education scenarios
human.mnlm.data <- subset(text.data, GPT == 0 & education == 0)
round(mean(human.mnlm.data$nchar))

# document-feature matrix for humans
human.dfm <- dfm_trim(dfm(tokens_ngrams(tokens(phrase(human.mnlm.data$mecabbed.text)), 1)), 
                      min_docfreq = 2)
human.dfm.matrix <- Matrix(as.matrix(human.dfm), sparse = TRUE)
dim(human.dfm.matrix)

## estimate the multinomial inverse regression model
GPT.mnlm.cov.statistics <- model.matrix(~ -1 + rating + gender + nationality + 
                                          statistics + disadvantaged + minimal + detailed, 
                                        GPT.mnlm.data)
GPT.mnlm.cov.needs <- model.matrix(~ -1 + rating + gender + nationality + 
                                     needs + disadvantaged + minimal + detailed, 
                                   GPT.mnlm.data)
human.mnlm.cov.statistics <- model.matrix(~ -1 + rating + gender + nationality + 
                                            statistics + disadvantaged + minimal + detailed, 
                                          human.mnlm.data)
human.mnlm.cov.needs <- model.matrix(~ -1 + rating + gender + nationality + 
                                       needs + disadvantaged + minimal + detailed, 
                                     human.mnlm.data)

cl <- makeCluster(16)
GPT.mnlm.results.statistics <- mnlm(cl, GPT.mnlm.cov.statistics, GPT.dfm.matrix)
GPT.mnlm.results.needs <- mnlm(cl, GPT.mnlm.cov.needs, GPT.dfm.matrix)
human.mnlm.results.statistics <- mnlm(cl, human.mnlm.cov.statistics, human.dfm.matrix)
human.mnlm.results.needs <- mnlm(cl, human.mnlm.cov.needs, human.dfm.matrix)
stopCluster(cl)

## identify words frequently used to justify decisions
# extract coefficients
GPT.mnlm.coef.statistics <- coef(GPT.mnlm.results.statistics)
GPT.mnlm.coef.needs <- coef(GPT.mnlm.results.needs)
human.mnlm.coef.statistics <- coef(human.mnlm.results.statistics)
human.mnlm.coef.needs <- coef(human.mnlm.results.needs)

# statistical scenarios in the ChatGPT survey
GPT.coef.rating.statistics <- GPT.mnlm.coef.statistics["rating", ]
GPT.specificity.statistics <- GPT.mnlm.coef.statistics["statistics", ]
GPT.freq.statistics <- colSums(GPT.dfm.matrix[GPT.mnlm.data$statistics == 1, ])
GPT.plot.coef.rating.statistics <- GPT.coef.rating.statistics[GPT.coef.rating.statistics != 0 & 
                                                                GPT.specificity.statistics > 0 & 
                                                                GPT.freq.statistics > 4]
GPT.plot.coef.rating.order.statistics <- order(abs(GPT.plot.coef.rating.statistics))
GPT.plot.coef.rating.statistics <- GPT.plot.coef.rating.statistics[GPT.plot.coef.rating.order.statistics]
GPT.plot.specificity.statistics <- GPT.specificity.statistics[GPT.coef.rating.statistics != 0 & 
                                                                GPT.specificity.statistics > 0 & 
                                                                GPT.freq.statistics > 4]
GPT.plot.specificity.statistics <- GPT.plot.specificity.statistics[GPT.plot.coef.rating.order.statistics]
GPT.plot.freq.statistics <- colSums(GPT.dfm.matrix[GPT.mnlm.data$statistics == 1, 
                                                   GPT.coef.rating.statistics != 0 & 
                                                     GPT.specificity.statistics > 0 & 
                                                     GPT.freq.statistics > 4])
GPT.plot.freq.statistics <- GPT.plot.freq.statistics[GPT.plot.coef.rating.order.statistics]

# customer needs scenarios in the ChatGPT survey
GPT.coef.rating.needs <- GPT.mnlm.coef.needs["rating", ]
GPT.specificity.needs <- GPT.mnlm.coef.needs["needs", ]
GPT.freq.needs <- colSums(GPT.dfm.matrix[GPT.mnlm.data$needs == 1, ])
GPT.plot.coef.rating.needs <- GPT.coef.rating.needs[GPT.coef.rating.needs != 0 & 
                                                      GPT.specificity.needs > 0 & 
                                                      GPT.freq.needs > 4]
GPT.plot.coef.rating.order.needs <- order(abs(GPT.plot.coef.rating.needs))
GPT.plot.coef.rating.needs <- GPT.plot.coef.rating.needs[GPT.plot.coef.rating.order.needs]
GPT.plot.specificity.needs <- GPT.specificity.needs[GPT.coef.rating.needs != 0 & 
                                                      GPT.specificity.needs > 0 & 
                                                      GPT.freq.needs > 4]
GPT.plot.specificity.needs <- GPT.plot.specificity.needs[GPT.plot.coef.rating.order.needs]
GPT.plot.freq.needs <- colSums(GPT.dfm.matrix[GPT.mnlm.data$needs == 1, 
                                              GPT.coef.rating.needs != 0 & 
                                                GPT.specificity.needs > 0 & 
                                                GPT.freq.needs > 4])
GPT.plot.freq.needs <- GPT.plot.freq.needs[GPT.plot.coef.rating.order.needs]

# statistical scenarios in the human survey
human.coef.rating.statistics <- human.mnlm.coef.statistics["rating", ]
human.specificity.statistics <- human.mnlm.coef.statistics["statistics", ]
human.freq.statistics <- colSums(human.dfm.matrix[human.mnlm.data$statistics == 1, ])
human.plot.coef.rating.statistics <- human.coef.rating.statistics[human.coef.rating.statistics != 0 & 
                                                                    human.specificity.statistics > 0 & 
                                                                    human.freq.statistics > 4]
human.plot.coef.rating.order.statistics <- order(abs(human.plot.coef.rating.statistics))
human.plot.coef.rating.statistics <- human.plot.coef.rating.statistics[human.plot.coef.rating.order.statistics]
human.plot.specificity.statistics <- human.specificity.statistics[human.coef.rating.statistics != 0 & 
                                                                    human.specificity.statistics > 0 & 
                                                                    human.freq.statistics > 4]
human.plot.specificity.statistics <- human.plot.specificity.statistics[human.plot.coef.rating.order.statistics]
human.plot.freq.statistics <- colSums(human.dfm.matrix[human.mnlm.data$statistics == 1, 
                                                       human.coef.rating.statistics != 0 & 
                                                         human.specificity.statistics > 0 & 
                                                         human.freq.statistics > 4])
human.plot.freq.statistics <- human.plot.freq.statistics[human.plot.coef.rating.order.statistics]

# customer needs scenarios in the human survey
human.coef.rating.needs <- human.mnlm.coef.needs["rating", ]
human.specificity.needs <- human.mnlm.coef.needs["needs", ]
human.freq.needs <- colSums(human.dfm.matrix[human.mnlm.data$needs == 1, ])
human.plot.coef.rating.needs <- human.coef.rating.needs[human.coef.rating.needs != 0 & 
                                                          human.specificity.needs > 0 & 
                                                          human.freq.needs > 4]
human.plot.coef.rating.order.needs <- order(abs(human.plot.coef.rating.needs))
human.plot.coef.rating.needs <- human.plot.coef.rating.needs[human.plot.coef.rating.order.needs]
human.plot.specificity.needs <- human.specificity.needs[human.coef.rating.needs != 0 & 
                                                          human.specificity.needs > 0 & 
                                                          human.freq.needs > 4]
human.plot.specificity.needs <- human.plot.specificity.needs[human.plot.coef.rating.order.needs]
human.plot.freq.needs <- colSums(human.dfm.matrix[human.mnlm.data$needs == 1, 
                                                  human.coef.rating.needs != 0 & 
                                                    human.specificity.needs > 0 & 
                                                    human.freq.needs > 4])
human.plot.freq.needs <- human.plot.freq.needs[human.plot.coef.rating.order.needs]

## Figures S7 and S8
png("Figure_S7.png", 
    width = 6, height = 2, unit = "in", pointsize = 6, res = 2000)
layout(matrix(1:2, 1, 2, byrow = TRUE))
par(mar = c(3.5, 3.5, 2, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = range(human.plot.coef.rating.statistics), 
     ylim = range(human.plot.specificity.statistics), 
     main = "Sphere and mechanism: Public (statistical)", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
text(human.plot.coef.rating.statistics, human.plot.specificity.statistics, 
     names(human.plot.specificity.statistics), 
     cex = sqrt(log10(human.plot.freq.statistics - 3)), 
     col = gray(ifelse(1 - abs(human.plot.coef.rating.statistics) < 0, 
                       0, 1 - abs(human.plot.coef.rating.statistics)), 0.5))
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Coefficient of ratings", 1, 2.5, font = 2)
mtext("Coefficient of the treatment condition", 2, 2.5, font = 2)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = range(human.plot.coef.rating.needs), 
     ylim = range(human.plot.specificity.needs), 
     main = "Sphere and mechanism: Public (customer needs)", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
text(human.plot.coef.rating.needs, human.plot.specificity.needs, 
     names(human.plot.specificity.needs), 
     cex = sqrt(log10(human.plot.freq.needs - 3)), 
     col = gray(ifelse(1 - abs(human.plot.coef.rating.needs) < 0, 
                       0, 1 - abs(human.plot.coef.rating.needs)), 0.5))
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Coefficient of ratings", 1, 2.5, font = 2)
mtext("Coefficient of the treatment condition", 2, 2.5, font = 2)
dev.off()

png("Figure_S8.png", 
    width = 6, height = 2, unit = "in", pointsize = 6, res = 2000)
layout(matrix(1:2, 1, 2, byrow = TRUE))
par(mar = c(3.5, 3.5, 2, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = range(GPT.plot.coef.rating.statistics), 
     ylim = range(GPT.plot.specificity.statistics), 
     main = "Sphere and mechanism: Public (statistical)", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
text(GPT.plot.coef.rating.statistics, GPT.plot.specificity.statistics, 
     names(GPT.plot.specificity.statistics), 
     cex = sqrt(log10(GPT.plot.freq.statistics - 3)), 
     col = gray(ifelse(1 - abs(GPT.plot.coef.rating.statistics) < 0, 
                       0, 1 - abs(GPT.plot.coef.rating.statistics)), 0.5))
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Coefficient of ratings", 1, 2.5, font = 2)
mtext("Coefficient of the treatment condition", 2, 2.5, font = 2)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = range(GPT.plot.coef.rating.needs), 
     ylim = range(GPT.plot.specificity.needs), 
     main = "Sphere and mechanism: Public (customer needs)", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
text(GPT.plot.coef.rating.needs, GPT.plot.specificity.needs, 
     names(GPT.plot.specificity.needs), 
     cex = sqrt(log10(GPT.plot.freq.needs - 3)), 
     col = gray(ifelse(1 - abs(GPT.plot.coef.rating.needs) < 0, 
                       0, 1 - abs(GPT.plot.coef.rating.needs)), 0.5))
axis(1, lwd = 0.5)
axis(2, lwd = 0.5)
mtext("Coefficient of ratings", 1, 2.5, font = 2)
mtext("Coefficient of the treatment condition", 2, 2.5, font = 2)
dev.off()

## extract representative texts for justification
# humans' justification (Table S5)
human.mnlm.data.statistics <- subset(human.mnlm.data, statistics == 1 & nchar > mean(nchar))
human.dfm.matrix.statistics <- human.dfm.matrix[human.mnlm.data$statistics == 1 & 
                                                  human.mnlm.data$nchar > mean(human.mnlm.data$nchar), ]
human.srproj.results.statistics <- srproj(human.mnlm.results.statistics, human.dfm.matrix.statistics)

head(human.mnlm.data.statistics$text[order((human.srproj.results.statistics[, "rating"] * 
                                              human.srproj.results.statistics[, "statistics"]))], 5)

human.mnlm.data.needs <- subset(human.mnlm.data, needs == 1 & nchar > mean(nchar))
human.dfm.matrix.needs <- human.dfm.matrix[human.mnlm.data$needs == 1 & 
                                             human.mnlm.data$nchar > mean(human.mnlm.data$nchar), ]
human.srproj.results.needs <- srproj(human.mnlm.results.needs, human.dfm.matrix.needs)

head(human.mnlm.data.needs$text[order((human.srproj.results.needs[, "rating"] * 
                                         human.srproj.results.needs[, "needs"]))], 5)

# ChatGPT's justification (Table S6)
GPT.mnlm.data.statistics <- subset(GPT.mnlm.data, statistics == 1)
GPT.dfm.matrix.statistics <- GPT.dfm.matrix[GPT.mnlm.data$statistics == 1, ]
GPT.srproj.results.statistics <- srproj(GPT.mnlm.results.statistics, GPT.dfm.matrix.statistics)

head(GPT.mnlm.data.statistics$text[order((GPT.srproj.results.statistics[, "rating"] * 
                                            GPT.srproj.results.statistics[, "statistics"]))], 1)

GPT.mnlm.data.needs <- subset(GPT.mnlm.data, needs == 1)
GPT.dfm.matrix.needs <- GPT.dfm.matrix[GPT.mnlm.data$needs == 1, ]
GPT.srproj.results.needs <- srproj(GPT.mnlm.results.needs, GPT.dfm.matrix.needs)

head(GPT.mnlm.data.needs$text[order(GPT.srproj.results.needs[, "rating"] * 
                                      GPT.srproj.results.needs[, "needs"])], 1)